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Abstract 

A new classification of sandpile models into universality classes is presented. 
On the basis of extensive numerical simulations, in which we measure an ex- 
tended set of exponents, the Manna two state model [S. S. Manna, J. Phys. A 
24, L363 (1991)] is found to belong to a universality class of random neighbor 
models which is distinct from the universality class of the original model of 
Bak, Tang and Wiesenfeld [P. Bak, C. Tang and K. Wiensenfeld, Phys. Rev. 
Lett. 59, 381 (1987)]. Directed models are found to belong to a universality 
class which includes the directed model introduced and solved by Dhar and 
Ramaswamy. 
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The introduction of sandpile models as a paradigm of self-organized criticality by Bak, 
Tang and Wiesenfeld (BTW) |I| stimulated numerous theoretical and numerical stud- 
ies 0-0. In these models, which are defined on a lattice, grains are deposited randomly 
until the height at some site exceeds a threshold, and becomes unstable. "Sand" is then 
distributed to the nearest neighbors. As a result of this relaxation process neighboring sites 
may become unstable, resulting in a cascade of relaxations called an avalanche. It was 
observed that these models are self-driven into a critical state which is characterized by a 
set of exponents fl|]. These include exponents that describe the distribution of quantities 
such as avalanche size and lifetime, and exponents which relate these properties of the dy- 
namics. Large scale simulations of the BTW model and some variants of it |||| were 
performed. The BTW model and the Manna two-state model were concluded to belong 
to the same universality class ||. Christensen and Olami later introduced an extended set 
of exponents 0. They measured the values of these exponents for the BTW model, and 
gave theoretical predictions and heuristic arguments for the values of some of the exponents. 
Continuous height models were also studied |IIJ and some aspects of universality were ex- 
amined A sandpile model with a preferred direction was introduced and solved by Dhar 



and Ramaswamy ]3[]. 

In this paper we present simulation results which suggest a new classification of sand- 
pile models into universality classes. The Manna two-state model is found to belong to a 
universality class of random relaxation models which is distinct from the BTW universality 
class. We first describe the different models, and define the properties of avalanches, with 
the exponents characterizing them. The models are defined onad dimensional lattice of 
linear size L. Each site is assigned a dynamic variable E(i) which represents some physical 
quantity such as energy, stress etc. In a critical height model a configuration {E(i)} is called 
stable if for all sites E(i) < E c , where E c is a threshold value. The evolution between stable 
configurations is by the following rules: 

(i) Adding energy. Given an arbitrary stable configuration {E(j)} we select a site i at 
random and increase E(i) by some amount 5E. When an unstable configuration is reached 



rule (ii) is invoked. 

(ii) The relaxation rule. If the dynamical variable at site i exceeds the threshold E c , 
relaxation takes place, whereby energy is distributed in the following way: 

£(i)^(i)-^A£(e) 

e 

(1) 

E(i + e) -f £(i + e) + A£(e), 

where e are a set of (unit) vectors from the site i to some neighbors. As a result of the 
relaxation the dynamic variable in one or more of the neighbors may exceed the threshold. 
The relaxation rule is then applied until a stable configuration is reached. The sequence of 
relaxations is an avalanche which propagates through the lattice. 



The parameters 5E and E c are irrelevant to the scaling behavior pJTT||. Thus the only 
factor determining the exponents is the vector AE, to be termed relaxation vector. For a 
square lattice with relaxation to nearest neighbors it is of the form AE = (E N , E E , E s , E w ), 
where E^ for example is the amount transferred to the northern nearest neighbor. The orig- 
inal BTW model is given by the vector (1, 1, 1, 1). The relaxation in the directed model of 
Dhar and Ramaswamy p3] is specified by any vector with ones in two adjacent directions 
and zeroes in the two other directions, such as (0,0,1,1). In a random relaxation model 
a set of neighbors is randomly chosen for relaxation. Such a model is specified by a set 
of relaxation vectors, each vector being assigned a probability for its application. As an 
example, a possible realization of a two-state model makes use of the six relaxation vectors 
(1,1,0,0), (1,0,1,0), (1,0,0,1), (0,1, 1,0), (0,1,0,1) and (0,0,1,1), each one applied with a proba- 
bility of 1/6. In Manna's two-state model || the variable is decreased to zero on relaxation, 
with sand distributed randomly among the nearest neighbors. We define a current 

J[A£]=5>£(e)e, (2) 

e 

which is the net flow in a relaxation. We also define 

3 = J2mE]P(AE), (3) 

AE 



which is the current averaged over the ensemble of relaxation vectors. Models can be clas- 
sified according to the value of the current J[Ai£] and its average, J. A model is called 
non- directed if J [AS] = (the BTW model for example). Random relaxation models such 
as the Manna two-state model, which satisfy J [AS] 7^ and J = 0, are called non-directed 
on average. Models with J 7^ are called directed. In this paper we present evidence that 
this is a classification into universality classes. 

Avalanches have various properties which can be measured in a simulation: size, area, 
lifetime, linear size, and perimeter. The size (s) of an avalanche is the total number of 
relaxation events that occured in the course of a single avalanche. The area (a) is the 
number of sites in the lattice where relaxation occured. Relaxation of all sites which exceed 
the threshold at a given time is considered a single time step. The lifetime (t) of an avalanche 
is the number of such steps. As for the linear size of an avalanche, there is no unique choice. 
A possible choice is the maximal distance (d) between the origin of the avalanche to sites of 
the avalanche cluster. Another possibility is the radius of gyration (r) of the cluster of sites 
where relaxation occured. A site belonging to the cluster of sites visited by an avalanche is 
defined to be a perimeter site if it has a nearest neighbor where no relaxation took place. The 
perimeter (p) is the number of perimeter sites. Thus we have a set of variables {s, a, t, r, d, p] 
which characterize an avalanche. The avalanche variables have probability functions which 
are assumed to fall off with a power law defined by P(x) ~ x l ~ Tx , where x G {s, a, t, r, d, p}. 
These variables also scale against each other in the form 

y~x lvx , (4) 

for x, y G {s,a,t,r,d,p}. The exact definition of the 7's is in terms of conditional expecta- 
tions values: S[y|x] ~ x lyx ||. The exponents are not independent. Scaling relations are 
found in J7[. We just note that 

7yx — Jxy , 

(5) 

'fzx = 'Jzy'fyx- 
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Avalanches are proven to be compact for BTW type models J?| but have a fractal boundary. 
It is reasonable to assume that the fractal dimension Df of the boundary is given by the 
scaling of the perimeter (p) against the linear size of the avalanche. It seems that for models 



which are non-directed the radius of gyration is the proper measure of size ||11|| . Therefore 
we identify Df with 7 pr . For directed models the maximum distance from the origin to the 
perimeter is the proper measure of size, and Df is identified with 7^. It is accepted that 
the dynamical exponent z of non-directed models should be identified with 7 tr ||1 1|| . In the 
case of directed models we identify the dynamical exponent with 7^. 

Having defined the models, we now describe the simulations. We used open boundary 
conditions and system sizes up to 512 2 , with 5 million grains dropped, in two dimensions; in 
three dimensions system sizes were up to 112 3 , with 20 million grains dropped. An algorithm 
due to Grassberger || was used. We ascertained the dynamics has reached the critical state 
by applying Dhar's "burning algorithm" ||, or by starting with a configuration belonging 
to the critical state. Manna's and our own simulation results for the BTW model indicate 
that the distribution exponents are system size dependent, with a logarithmic convergence 
to the infinite system values. The values of the 7's on the other hand, seem to be almost 
independent of system size. Moreover, we found that the relations that specify the 7's hold 
during avalanches as well, and are not just a scaling property of completed avalanches. Thus 
the 7's provide a robust characterization of the dynamical properties of a sandpile model, 
and can be used for a reliable classification of sandpile models into universality classes. 

Previous studies clearly show that directed and non-directed models belong to different 
universality classes On the basis of Manna's simulation results it was concluded 

that the Manna two-state model and the BTW model are in the same universality class 
||. This conclusion is based on measurements of a limited set of exponents: T s ,r t and j ts . 
We measured the extended set of exponents introduced by Christensen and Olami, and the 
fractal dimension. The 7's we obtained in two dimensions are listed in Table |. Our results 
are consistent with known analytical results and simulation data: Dhar and Ramaswamy's 
analytical solution of a directed model [||]; simulation results and scaling arguments given by 



Christensen and Olami simulation results of Manna A momentum-space analysis 

of a Langevin equation indicates that for the BTW model z = (2 + d)/3 [[II]. Our results 
for j r t which is identified with 1/z, confirm this scaling relation. This agreement supports 
our observation that the 7's are size independent, and indicates that we are in the right 
avalanche size regime for the observation of 7 rt . On the basis of the difference in the 7's 
for the BTW and two-state models we conclude that the two models are not in the same 
universality class (Fig. [l]). 

In order to establish that the classification introduced above is a classification into uni- 
versality classes we provide evidence that some details of the models are irrelevant (Fig. [|). 
Simulation results of the BTW model on the triangular lattice and square lattice were com- 
pared @,[n] . No significant difference was reported. We define N as the number of states of 
the E(i) in stable configurations of discrete models. When the components of the relaxation 
vector are all l's, N also equals the number of neighbors. In sandpile models the question 
of the lattice dependence or interaction range dependence of the exponents is actually a 
question of the dependence on N. We observed a crossover effect when increasing N. The 
scaling obtained for the BTW model on a square lattice (N=4) is shifted to larger avalanches 
when N is increased. Similar cross-over was observed in the other universality classes. Note 
that the requirement that J [AS] = does not imply isotropy. This is the reason the uni- 
versality class was called non- directed, rather than isotropic. As an example, a model with 
a toppling vector (1,2,1,2) fulfills this requirement, and simulations show that it belongs to 
the universality class of non-directed models. 

Continuous models were simulated as well. There are two types of realizations of con- 
tinuous models. In one, the variables are turned into continuous variables, and when the 
amount of sand added is not a multiple of the amount distributed on relaxation (or is a 
random variable taking such values) then the height profile is turned into a continuous dis- 
tribution. The other is the Zhang realization, where on relaxation the dynamic variable 
is decreased to zero and sand distributed equally among the nearest-neighbors [|I(J. Both 
types seem to be in the same universality class |yj . This is indicated by our simulations as 
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well. 

There is a number of possible realizations of a two-state model. The neighbors to which 
sand is distributed can be chosen as distinct (no neighbor chosen twice) or not. In Manna's 
two-state model || the variable is decreased to zero on relaxation, with sand distributed 
randomly among the nearest neighbors. In this case the relaxation process depends on the 
variable value. Continuous variants of the model may also be defined. We have simulated 
realizations of such models and all were in the same universality class (Fig. 0). Simulations 
of two-state models were performed with annealed randomness only [|P2jl . 



On the basis of the wave structure of avalanches in the BTW model ||13|| , it can be shown 
that avalanches have a "shell" structure, i.e. the sites which relaxed at least n + 1 times form 
a connected cluster with no holes which is contained in the cluster of sites which relaxed 
at least n times (Fig. 0(a)). Avalanches in random relaxation models do not share this 
property, and their structure is more irregular. A typical avalanche in a two-state model 
is shown in Fig. [|(b) . These geometrical differences reflect in the fractal dimension of the 
boundary, which is greater for the two-state model. 

The distinction between the universality classes of non-directed models and models which 
are non-directed only on average holds in three dimensions as well (Table ||). The difference 
is less marked because the exponents are nearing their mean field values. 

Directed models also form a universality class. In addition to the models studied by 
Dhar and Ramaswamy, where the relaxation vector is of the form (1,1,0,0) or (1,1,1,0), we 
simulated models with the relaxation vectors (1,1,1,2) and (1,1,2,2). In the latter, multiple 
relaxations are possible, but it does not reflect in the scaling behavior. We found the same 
exponent values in all these models. The values we obtained in simulations (Table |) are 
in agreement with the analytical solution. Directed models with a random relaxation rule 
show cross-over. 

In summary, using extensive numerical simulations we identified three universality classes 
in sandpile models: (a) non-directed models (BTW model); (b) random relaxation models 
which are non-directed only on average (Manna two-state model) and (c) directed models 



(Dhar and Ramaswamy's directed model). These universality classes correspond to a classifi- 
cation according to the value of the current J[A.E] and its average. Locally non- conservative 
models also show SOC and may form a different universality class M. 



Recently Pietronero et al. jnj introduced a novel theoretical framework for calculating 
the exponents of sandpile models, in a manner which immediately reveals their universality. 
Within their scheme, which is purely phenomenological, the Manna two-state model and 
the BTW model are found to be in the same universality class. Its failure to distinguish 
between the two models indicates that some key ingredient is missing from their scheme. We 
suspect that multiple relaxation is the missing element. Work is now in progress to extend 
the procedure to include some form of multiple relaxation. 

Acknowledgements We thank Z. Olami, P. Bak, E. Domany, G. Grinstein, C. 
Jayaprakash, T. Kaplan and M. Paczuski for helpful discussions. 
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FIGURES 

FIG. 1. Simulation results for the BTW model (circles) and two-state model (squares). -E[r|t] 
(average radius of gyration for given avalanche lifetime) vs. t is displayed in (a), yielding j r t. In 
(b) we show a graph of i?[s[a] vs. a which yields 7 sa - Their values are listed in Table |. Data was 
binned with bin size increasing exponentially. System size is 512 2 , with 10 7 grains dropped. These 
results indicate that the two models belong to different universality classes. 

FIG. 2. Simulation results showing the universality of models in the BTW (non-directed) and 
two-state (non-directed on average) universality classes. Graph shows 2£[s|a] vs. a for system size 
128 2 . Unless otherwise stated simulations were performed on a square lattice. For models in the 
BTW class we see data collapse on a curve with ^ sa = 1.06. For the random relaxation models 
7sa = 1.23 ±0.01. 

FIG. 3. Typical avalanche structure for the BTW model (a) and two-state model (b). 
Grey-scale indicates the number of relaxations which occurred at each site during an avalanche. 
White represents zero relaxations, and black represents the maximal no. of relaxations (10 in (a), 
45 in (b)). System size is 150 2 . Note the shell structure in the BTW avalanche (an analytically 
provable property) vs. the irregular structure of the avalanche in the two-state model. These qual- 
itative geometrical differences translate into quantitative differences in exponent value, especially 
the fractal dimension of the boundary. 
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TABLES 

TABLE I. 7 exponents for universality classes in two dimensions. The other 7's can be found 
from the scaling relations, Eq. ^ . The values of these exponents were observed to be independent of 
system size. The typical spread of data for different runs of different models within the universality 
class is ±0.01 about the mean. 

TABLE II. Exponents in three dimensions for the BTW model and a three-state random 
relaxation model (non-directed on average). Distribution exponents are given for system size of 
96 3 . 
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Exponent 




model 






BTW 


two-state 


Directed 




0.76 


0.67 


1.00 


1st 


1.62 


1.70 


1.51 


lot 


1.53 


1.35 


1.51 


Isa 


1.06 


1.23 


1.00 


Df 


1.26 


1.42 


1.00 


a In non-directed models z 


is identified with jtr, 


and in directed models it 


is identified with 7^. 
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Exponent BTW 3-state 

t s 2.35 2.43 

r a 2.35 2.46 

7rt (1/z) 0.60 0.54 

lst 1.78 1.80 

lat 1.78 1.72 

Isa 1.00 1.06 
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